
le = [73.35
73.7
74.25
74.55
74.8
75.05
75.5
75.85
75.95
76.15
76.35
76.5
76.55
76.75
76.9
77
77.1
77.25
77.25
77.4
77.6
77.75
77.85
78
78.15
78.25
78.5
78.65
78.85
79.05
79.25
79.4
79.55
79.7
79.9
80.05
80.2
80.4
80.5
80.7
80.85
80.95
81.1
81.2
81.3
81.4
81.55
81.65
81.75
81.85
81.95
82.15
82.25
82.35
82.5
82.6
82.7
82.8
82.85
82.95
83.05
83.15
83.25
83.3
83.4
83.5
83.6
83.65
83.75
83.85
83.95
84
84.1
84.2
84.3
84.3
84.4
84.5
84.6
84.65
84.75
84.85
84.9
84.95
85.05
85.15
85.2
85.3
85.4
85.45
85.5
85.6
85.65
85.75
85.85
85.85
85.95
86.05
86.1
86.2
86.2
86.3
86.4
86.45
86.5
86.6
86.65
86.75
86.75
86.85
86.95
87
87.05
87.1
87.2
87.2
87.3
87.4
87.4
87.5
87.55
87.6
87.65
87.75
87.75
87.85
87.95
87.95
88.05
88.1
88.15
88.2
88.3
88.3
88.4
88.45
88.5
88.55
88.6
88.65
88.75
88.75
88.85
88.85
88.95
89
89.05
89.1
89.15
89.2
89.25] ;

le_year = 1940:2090 ;

%%

kappa0 = 2.2 ; %1
kappa1 = 15 ; %15
kappa2 = 0.68; %1.25
kappa3 = 0.03 ; %0.25

le_1 = le(le_year==str2num(year_)) ;

fun = @(x) (1/(sqrt(2*pi)*kappa3*le_1)) .* exp(-(1/2) .* ((x-kappa2*le_1)./(kappa3*le_1)).^2) ;

for ii=1:80
    
    vvec(ii) = (kappa1*ii / le_1) ;
    
    vvec(ii) = kappa0 + vvec(ii) * integral(fun,-Inf,ii) ;
    
end

%hold on; plot(16:95,vvec,'r')

%%

tvv = 0 ;

le_vec = le ; %(1-tvv)*le_1 + tvv*
le_vec = (1-tvv)*le_1 + tvv*le ;

kappa0_vec = (1-tvv)*kappa0+tvv*[linspace(kappa0,2,74) linspace(2, 2, 131-74)];%2
kappa1_vec = (1-tvv)*kappa1+tvv*[linspace(kappa1,40,74) linspace(40, 40, 131-74)];%150
kappa2_vec = (1-tvv)*kappa2+tvv*[linspace(kappa2,1,74) linspace(1, 1, 131-74)];%linspace(kappa2,0.65,131) ;%150
kappa3_vec = (1-tvv)*kappa3+tvv*[linspace(kappa3,0.3,74) linspace(0.3, 0.3, 131-74)];%linspace(kappa3,0.35,131) ;

kappa4_vec = zeros(80,131) ;
kappa4_vec(1,:) = (1-tvv)*0 + tvv*linspace(3.0,6.0,131) ;

for jj=1:131
    for ii=2:80
        kappa4_vec(ii,jj) = kappa4_vec(ii-1,jj)*0.65 ;
    end
    
    %kappa4_vec(9,:) = linspace(-0.5,-0.5,131) ;
    %kappa4_vec(10,:) = linspace(-0.6,-0.6,131) ;
    %kappa4_vec(11,:) = linspace(-0.7,-0.7,131) ;
    %kappa4_vec(12,:) = linspace(-0.8,-0.8,131) ;
    %kappa4_vec(13,:) = linspace(-0.9,-0.9,131) ;
    %kappa4_vec(14,:) = linspace(-1.0,-1.0,131) ;
    %kappa4_vec(15,:) = linspace(-1.1,-1.1,131) ;
    
    %for ii=16:80
    %    kappa4_vec(ii,jj) = kappa4_vec(ii-1,jj)*0.9 ;
    %end    
end

for jj=1:131

fun = @(x) (1/(sqrt(2*pi)*kappa3_vec(jj)*(le_vec(jj)))) ...
    .* exp(-(1/2) .* ((x-kappa2_vec(jj)*(le_vec(jj)))./(kappa3_vec(jj)*(le_vec(jj)))).^2) ;

for ii=1:80
    
    vvec_sh(ii,jj) = (kappa1_vec(jj)*ii / (le_vec(jj))) ;
    
    vvec_sh(ii,jj) = kappa4_vec(ii,jj) + kappa0_vec(jj) + vvec_sh(ii,jj) * integral(fun,-Inf,ii) ;
    
end

end

vvec_end = vvec_sh ; %vvec_end = repmat(vvec',[1 131]) ;

tmp = linspace(0,0,131) ;

vvec_sh = 1*vvec_sh - 1*repmat(vvec',[1 131]) + repmat(tmp, [80,1]) ;